Section: New Results

Development of numerical methods

Participants : Aurore Back, Nicolas Besse, Jean-Philippe Braeunig, Anaïs Crestetto, Emmanuel Frénod, Philippe Helluy, Sever Hirstoaga, Ahmed Ratnani, Eric Sonnendrücker.

Application of isogeometric analysis to plasma physics simulations

Mainly around the PhD thesis of Ahmed Ratnani [13] which has been defended in October 2011, we have been using the concept of isogeometric analysis introduced by Hughes and co-workers [77] which consists in representing the computational domain as well as the numerical solution of the equations with NURBS (Non Uniform Rational B-Splines).

In [26] we introduced a time-domain conforming Finite Element solver using arbitrary order B-Splines as basis functions. The discrete function spaces used in the Finite Element formulation form a De Rham sequence, which has proved to be an important property for numerical Maxwell solvers. In particular, they allow to have a simple relation between spline coefficients of the magnetic and electric field that is independent on order and geometry for one of Ampère's or Faraday's law. The other then necessarily involves a discrete Hodge operator which depends on order and geometry. High-order energy conserving leap-frog schemes have been validated with this solver.

In [21] , we developed an arbitrary order B-Spline Finite Element solver for the quasi-neutrality equation that is generally coupled to gyrokinetic Vlasov-solvers. Compared to the previous solver used in GYSELA which was spectral in the angular variable and second order Finite Differences in the radial variable this solver can be of high-order in both direction. This enables us for a given accuracy to decrease the number of grid points. Moreover, thanks to the periodicity in the angular variable and the tensor product structure of the problem, we could introduce a fast diagonalization method using a FFT such that the cost of the new solver only marginally depends on the order and is only slightly higher than the cost of the previously used method. Another important new algorithm introduced in this work is the decoupling of the parallel and transverse parts of the equation, by solving successively for the average value along the parallel direction and the remaining part.

In [28] we present an axisymmetric PIC code based on isogeometric analysis, which was initially the IsoPIC project (supported by CEA Gramat) of CEMRACS 2010. The goal of this study is to use it for solving the system of Vlasov-Maxwell equations. The idea is to develop an axisymmetric Finite Element PIC (Particle-In-Cell) code in which specific spline Finite Elements are used to solve the Maxwell equations (in 2D transverse electric mode) and the same spline functions serve as shape function for the particles. The computational domain itself is defined using splines or NURBS. We are in particular interested in the emission of electrons in a diode with hemispherical cathode (thanks to symmetry in $\theta$ direction, we can consider the 2D axisymmetric geometry).

Spline discrete differential forms

In [11] , [36] we have developed the concept of spline discrete differential forms which can be used to discretized equations defined using the notions of exterior calculus. These have been applied for the numerical solution of the Maxwell and the Vlasov equations. Hodge operators either using a dual grid or a weak formulation are developed and commuting diagram properties are proved.

Drift-kinetic simulations

We implemented the conservative semi-Lagrangian method in the GYSELA code which is based on the classical backward semi-Lagrangian method which is not exactly conservative. We noticed that for the conservative method it is essential that the advection field remains numerically exactly divergence free in order to avoid numerical instabilities. In addition specific limiters for the conservative method were developed and comparison between the backward semi-Lagrangian method, the conservative semi-Lagrangian method with 1D splitting and the same method with an unsplit Finite Volume like formulation in the $(r,\theta )$ plane which provides a better conservation of volume [43] .

Waterbag simulations

In [31] we apply the multi-water-bag model and the method of moments to the Vlasov-Poisson system in a case where the solution becomes multivalued. The motivation of this study is that the kinetic Vlasov-Poisson model is very expensive to solve numerically. It can be approximated by a multi-water-bag model in order to reduce the complexity. This model amounts to solve a set of Burgers equations, which can be done easily by finite volume methods. However, the physical solution can become multivalued (filamentation appears). In this case, shocks appear in the simulation and we lose information about the filaments. To catch them, we can use a moment method. We describe here the two models and present several numerical experiments.

A linear analysis code CYLGYR based on the gyrowaterbag model in cylindrical geometry has been developed. It enables, starting from a given equilibrium configuration, to obtain the whole set of modes that can exist. It was used to validate the linear phase of the previously developed non linear semi-lagrangian gyrowaterbag code in cylindrical geometry GMWB3D-SL [61] . Excellent agreement for the growth rates (eigenvalues) and the radial envelopes (eigenfunctions) has been obtained for global eigenmodes. On the other hand the linear code CYLGYR gives results in excellent agreement with those given by the linear kinetic code KINEZERO. We are now using as well the linear gyrowaterbag code CYLGYR and the non linear code GMWB3D-SL to compare the quasi-linear and non-linear fluxes and thus measure the validity of the quasi-linear approach for gyrokinetic turbulence. Such comparisons will also be performed with the gyrokinetic code GYSELA in cylindrical geometry.

Validation of the quasi-linear theory

We have developed and optimized a parallel semi-Lagrangian code for the numerical resolution of the Hamiltonian Vlasov-wave model in two phase-space dimensions. Using this code to perform a statistical study on a large number of runs of the system we have showed that the quasi-linear theory, whose aim is to justify the approximation of a self-coherent hamiltonian system like the Vlasov-Poisson model by diffusive self-coherent model of Fokker-Planck type, was valid in the strongly chaotic non linear regime of 1D electrostatic turbulence provided it is regarded from a statistical point of view [17] .